/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/***************************************************************************************************************************************
 *
 * BRAHMS - a Biomembrane Reduced-ApproacH Molecular Simulator
 *   
 * Brahms is a molecular dynamics program specifically designed for the simulation of lipid and water systems modeled with the ELBA 
 * coarse-grain force field [Orsi & Essex, PLoS ONE, 6, e28637 (2011)].
 *
 * The development of Brahms has been funded by: 
 *  - UK taxpayers, through government funding agencies: BBSRC, EPSRC;
 *  - private companies: Johnson & Johnson, Unilever.
 *
 * The main author of Brahms is Mario Orsi (www.soton.ac.uk/~orsi).
 *
 * Many crucial components of Brahms, such as the main functions and data structures, the cell-subdivision/neighbour-list algorithm, 
 * and several useful macro definitions, have been developed following the excellent book 'The Art of Molecular Dynamics Simulation'
 * by D.C. Rapaport (2004, 2nd ed). Specific references to this book throughout the code are identified by the format [Rapaport, pp.XX-XX]. 
 *
 * Brahms is intended to be clear code, easy to read, develop and maintain. To achieve these objectives, 
 * a necessary condition involves consistent adherence to a set of style rules: 
 *  - as it is very common practice, constants should be uppercase words, possibly separated by underscores (THIS_IS_A_CONSTANT)
 *  - variables are "initial caps", except for the very first character: thisIsAVar
 *  - functions are "initial caps", including the very first character: ThisIsAFunc
 *  - macro definitions are also "initial caps": #define VScaleTwo(v, s) (v).x *= s, (v).y *= s
 *  - types are also "initial caps": typedef struct { real x, y, z; } VecR; 
 *  - pointer variables should have the "Ptr" suffix (*filePtr) unless the pointers are really arrays (like the site[] array)
 *  - all variables declarations should have a comment explaining what they do
 *  - variable names should be one ('variable'), two ('anotherVariable'), or more words ('yetAnotherUsefulVariable'). 
 *    In no cases should a single or double character variable name be allowed, unless the full name is one or two characters such as "pi."
 *  - the common indentation size should be 2 spaces(*)
 *  - when the indent causes the program to run into the right side of the screen (for example because of too many nested "ifs"), 
 *    consider breaking the code into simpler functions
 *  - curly braces should be on the same line as the statements
 *  - all functions should be preceded by a comment block explaining them. This comment block should stand out;
 *  - functions should be at most two or three pages long. Anything longer should probably be split up into smaller, simpler functions;
 *  - each case statement in a switch should end with "break" or the comment "fall through" (to let people know we didn't forget a break
 *  - reliability and readability are major targets - slow and steady beats complex, wonderful, and fast every time; 
 *    BUT: consider exceptions to this rule where speed is crucial (as in the calculation of nonbonded interactions)
 *  - Emacs users: insert extra parentheses so that Emacs will indent the code properly. 
 *    For example, the following indentation looks nice if you do it by hand,
 *      v = rup->ru_utime.tv_sec*1000 + rup->ru_utime.tv_usec/1000
 *          + rup->ru_stime.tv_sec*1000 + rup->ru_stime.tv_usec/1000;
 *    but Emacs would alter it. Adding a set of parentheses produces something that looks equally nice, and which Emacs will preserve:
 *      v = (rup->ru_utime.tv_sec*1000 + rup->ru_utime.tv_usec/1000
 *           + rup->ru_stime.tv_sec*1000 + rup->ru_stime.tv_usec/1000);
 *   - **The most important rule**: If any of the above rules keep you from making a clear, easily to understand, readable, 
 *    reliable program, discard them. Be prepared to defend your decision, but don't let the rules get in the way of doing what's right.
 * Most of these rules are taken from Rapaport [pp.20-25], Oualline [Practical C Programming, O'Reilly, 3rd edition, p.38] and the GNU 
 * coding standards (http://www.gnu.org/prep/standards).
 *
 * Brahms manual, examples, and past and present code versions, can be found online: www.soton.ac.uk/~orsi/brahms 
 *
 * Scientific papers based on Brahms simulations (please cite these papers if you use or reference Brahms!):
 *  - Orsi, M. and J. W. Essex, The ELBA force field for coarse-grain modeling of lipid membranes, PLoS ONE, 6, e28637 (2011).
 *  - Orsi, M., M. G. Noro and J. W. Essex. Dual-resolution molecular dynamics simulation of antimicrobials in biomembranes, 
 *    Journal of the Royal Society Interface, 8, 826 (2011).
 *  - Orsi, M. and J. W. Essex. Permeability of drugs and hormones through a lipid bilayer: insights from dual-resolution 
 *    molecular dynamics, Soft Matter,  6, 3797 (2010). 
 *  - Orsi, M., J. Michel and J. W. Essex. Coarse-grain modelling of DMPC and DOPC lipid bilayers, Journal of Physics: 
 *    Condensed Matter, 22, 155106 (2010). 
 *  - Orsi, M., W. E. Sanderson and J. W. Essex. Permeability of small molecules through a lipid bilayer: a multiscale simulation study, 
 *    Journal of Physical Chemistry B, 113, 12019 (2009).
 *  - Orsi, M., D. Y. Haubertin, W. E. Sanderson and J. W. Essex. A quantitative coarse-grain model for lipid bilayers, 
 *    Journal of Physical Chemistry B, 112, 802 (2008). 
 *
 ****************************************************************************************************************************************/
 
/*****************************************************************************************************************************************
 * main.c -- this module contains the definition of the main global variables, together with functions and function-calls for setting up *
 *           the system, running the main molecular dynamics loop, computing various properties and outputting data of interest          *
 *****************************************************************************************************************************************/

#include "dataStructs.h"
#include "in_namelist.h"

/* external data definitions - global variables */
Site *site; /* array of Site structures, contains several site features (see definition in dataStructs.h) */
VecR *siteForce; /* site forces */
VecR *siteTorque; /* torques in the space-fixed frame */
VecR *siteCoords; // site coordinates (x,y,z)
VecR *siteVelocity; // site velocity
VecR *siteAngMom;  /* angular momentum of rigid-body sites in the body-fixed frame */
VecR *siteOrientVec; /* unit vector associated with rigid-body sites. Two types of rigid-body sites currently occur in Brahms:
			i) SSD water sites: nonsymmetric rigid bodies. The orientation vector corresponds to the electrical dipole vector 
			ii) Dipolar lipid sites (eg, glycerol and ester): symmetric rigid bodies. The orientation vector is oriented from 
			the bottom (tails) to the top (headgrouops) of the lipid molecule. Direction coincides with the embedded electrical
			dipole. It follows that a positive electrical dipole is aligned with "siteOrientVec", whereas a negative
			dipole is antialigned with it. */
RMat *siteRotMat; /* rotation matrix for rigid-body sites */
VecR *atomCoordinates; 
Atom *atom;           /* array of Atom, represents z-constrained solute */
VecR region;          /* edge lengths of the simulation region */ 
VecI initUcell;       // size of unit cell array for initial state - only used when starting new run
real rCutMax; // longest cutoff radius
real timeNow; /* current time */

// Properties of interest [Rapaport, pp.30-32]- REMEMBER to initialise Prop variables (AccumProps function)
Prop totEnergy; 
Prop potEnergy; 
Prop systemTemperature;
Prop lipTemperature;
Prop watTemperature;
Prop solTemperature;
Prop rotWatTemp; 
Prop transWatTemp; 
Prop rotLipTemp; 
Prop transLipTemp;
Prop pressure; /* P and average */
Prop boxVolume; 
Prop xyArea; // xy-plane
Prop lipVolume;
Prop lipArea;
Prop surfaceTension; // surface tension
Prop headGroupDipole; 
Prop zConstraintForce;
Prop dens;

int moreCycles; /* logical */
int nTypes; /* #site types, e.g., 5 ( water:0, chol:1, phos:2, gly:3, tail:4 )*/
int nDOPCsDSPCs; // #DOPC or DSPC lipids
int nDOPEs; // #DOPE lipids
int nLipids; /* #lipid molecules */
int nWaters; /* #water molecules - INPUT PARAMETER */
int nSites; /* #sites */ 
int nAtoms; // #atoms of the z-constrained solute
int nSolutes; // #solutes
int nPointMasses; /* #Lennard-Jones sites */
int nSymmRigBodies; /* #symmetric rigid bodies (e.g., dipolar sites) */
int nNonSymmRigBodies; // #non-symmetric rigid-bodies (e.g., SSD sites)
int stepCount; /* timestep counter */
int stepLimit; /* total run length */

real totalMass;
RMat boxMatrix;

/* cell subdivision & neighbour list: for efficient calculation of nonbonded pair interactions  */
real rNebrShell; // thickness of neighbor list shell
int nebrTabFac; // used to determine nebrTabMax [Rapaport, p. 513]
int nebrTabMax; // maximum neighbor-list length [Rapaport, p. 513]
int maxEdgeCells; // maximum size of cell array [Rapaport, p. 512]

/* lateral diffusion variables */
VecRTwo *xyLipidCom, *xyLipidComTrue; // array of 2D xy-vector for the mass centre of each lipid (and possibly solute)
VecRTwo xySoluteCom, xySoluteComTrue;

VecR *lipUnwrappedCoords; // array of "true", unwrapped, lip coords - for visualisation and correct lat COM computation 

/* wat rdf variable */
real *histWatRdf, rangeWatRdf;  
int countWatRdf, limitWatRdf, sizeHistWatRdf;

// local (lateral really) pressure profile variables
real *totalHistLpp, *histLppTau1, *histLppTau2, *histPress_xx, *histPress_yy, *histPress_zz; 
// *watHistLpp, *hgHistLpp, *glyHistLpp, *tailHistLpp, *vdwHistLpp, *electroHistLpp, *bondHistLpp, *dummyHistLpp;

// water polarisation profile variables
real *histWpp;
int *nWatsPerWppSlab;

// miscellany
int doCheckpoint; // logical - controls checkpointing [Rapaport, pp. 500-504]
int stepCheckpoint; // checkpointing frequency (number of timesteps between generations of checkpoint files)
int stepPdb; // output frequency for "trajectory.pdb"
int reCenterBilayer; 
int writeAreaVol;

int zConstraint; // switch: if set to 1 it turns on the z-contraint algorithm 
FILE *fPtrZCF; // file that records the zConstraintForce, it opened in setUpJob and closed at the end of main()
FILE *fPtrAV; // file that records the zConstraintForce, it opened in setUpJob and closed at the end of main()

/* Topology input parameters - if not supplied in Brahms units, remember to do the conversion */
real sigWat, epsWat, sigChol, sigAmine, sigPhos, sigGly, sigEst, epsChol, epsAmine, epsPhos, epsGly, sigCharge, epsCharge; 
real massWat, massChol, massAmine, massPhos, massGly, massEst, massTail;
real inertiaWat, inertiaEst, inertiaGly;
real epsTail, epsEst, estDipole, sigTail;
real cholCharge, phosCharge, glyDipole, watDipole;
real **sig, **eps, **sigSquared, *mass, *inertia;
real *rLipid; /* vector containing lipids Centre Of Geometry position */
  
extern VecI cells; 
extern real rCutWatWat, rCutSolute, rCutSoluteElse, rCutLipLip;
extern const int applyBarostat, flexBox;
extern const int diffusion, stepDiffuse, latDiff, stepLatDiff, writeLipLatMotion;
extern const int rdfWat, stepWatRdf, edp, stepEdp, sizeHistEdp, lpp, stepLpp;
extern const int sizeHistLpp, wpp, stepWpp, sizeHistWpp;
extern const int epp, stepEpp, sizeHistEpp;  
extern const real tauP, tauT, deltaT;
extern real lppDeltaZ, edpDeltaZ, wppDeltaZ, eppDeltaZ;

/* external function declarations: */
void CheckBaroThermoParameters( real, real, real );
void CheckMinimumImageConvention( real );
void CheckStepParameters( int );
void ComputeAndUnwrapSoluteLatCom(); 
void ComputeLipLatCom();
void ConvertInputParametersToBrahmsUnits(); 
void DumpLammpsData();
void DumpLammpsForcefield();
void EvalDiffusion();
void EvalEdp();
void EvalLatDiff();
void EvalLpp();
void EvalWpp( int, real );
void EvalEpp();
void EvalWatRdf();
void OutputFinalConfiguration();
void WriteLipLatMotion(); 
void PrintBrahmsVersion();
void PrintSoluteLatMotion();
void PutCheckpoint();
void PutConfig();
void SetLipidTopologies();
void SetupJob( FILE* );
void SingleStep();
void PrintNameLists( FILE* );
void PrintPdb( FILE* );
void GetNameLists( char** );
void SetupInterrupt();
void UnwrapLipLatCom(); 
void WriteAreaVolume();
void WriteLammpsData( FILE* );
void WriteLammpsForcefield( FILE* );

/***********************************************************************
 * SetParams -- sets various parameters; adapted from [Rapaport, p.31] *
 ***********************************************************************/
static void SetParams( FILE *fp ) 
{
  const real WATER_NUMBER_DENSITY = 33.46; /* [molecules/nm^3] - experimental number density of H20 */

  SetLipidTopologies(); 

  if ( !nLipids && ( nWaters>255 ) ) { // build "pure" system of water from scratch
    VSCopy( region, 1. / pow( WATER_NUMBER_DENSITY / 4., 1./3. ), initUcell );
    nSites = 4 * VProd( initUcell );
  } 
  
  ConvertInputParametersToBrahmsUnits( fp ); 

  boxMatrix.u[0]=region.x; boxMatrix.u[4]=region.y; boxMatrix.u[8]=region.z;
  boxMatrix.u[1]=boxMatrix.u[2]=boxMatrix.u[3]=boxMatrix.u[5]=boxMatrix.u[6]=boxMatrix.u[7]=0.;

  boxVolume.val = MatDet( boxMatrix.u );
  xyArea.val = region.x * region.y;

  rCutMax = Max( rCutSolute, Max3( rCutSoluteElse, rCutLipLip,  Max( rCutLipLip, rCutWatWat ) ) ); 
  CheckMinimumImageConvention( rCutMax );

  VSCopy( cells, 1. / ( rCutMax + rNebrShell ), region ); 
  maxEdgeCells = 1.3 * Max3( cells.x, cells.y, cells.z );
  nebrTabMax = nebrTabFac * nSites;  

  // checking input parameters
  CheckBaroThermoParameters( tauP, tauT, deltaT );
  CheckStepParameters( nLipids );
}

/************************************************************************************************************
 * main -- controls system's initialisation, evolution in time and analysis - adapted from [Rapaport, p.21] *
 ************************************************************************************************************/
int main( int argc, char **argv ) 
{ 
  FILE *fPtr;

  PrintBrahmsVersion();
  GetNameLists( argv );
  PrintNameLists( stdout );
  fPtr = fopen( "inputParameters.log", "w" ); 
  PrintNameLists( fPtr ); 
  fclose( fPtr );

  fPtr = fopen( "simulationSetup.log", "w" ); 
  SetParams( fPtr ); 
  SetupJob( fPtr ); 
  fclose( fPtr );

  DumpLammpsForcefield( fPtr ); 
  DumpLammpsData( fPtr );

  SetupInterrupt();
  moreCycles = 1;
  
  //  PutConfig();
  while ( moreCycles ) {
    SingleStep();
    if ( doCheckpoint && stepCount % stepCheckpoint == 0 ) PutCheckpoint();  
    if ( applyBarostat && flexBox && writeAreaVol ) WriteAreaVolume(); 
    if ( stepCount % stepPdb == 0 ) {
      fPtr = fopen( "trajectory.pdb", "a" ); 
      PrintPdb( fPtr ); 
      fclose( fPtr );
    }

    /* analyze various properties, output average values: */
    if ( latDiff && ( ( stepCount % stepLatDiff ) == 0 ) ) {
      ComputeLipLatCom();
      EvalLatDiff();
      if ( writeLipLatMotion ) {
	UnwrapLipLatCom(); 
	WriteLipLatMotion(); 
      }
    }
    if ( diffusion && ( ( stepCount % stepDiffuse ) == 0 ) ) EvalDiffusion();
    if ( rdfWat && ( ( stepCount % stepWatRdf ) == 0 ) ) EvalWatRdf();
    if ( edp ) {
      edpDeltaZ = region.z / sizeHistEdp; // computing the binSlab size (slab thickness) - used in edp.c
      if ( ( stepCount % stepEdp ) == 0 ) EvalEdp();
    }
    if ( lpp ) {
      lppDeltaZ = region.z / sizeHistLpp; // computing the binSlab size (slab thickness) - used in lpp.c
      if ( ( stepCount % stepLpp ) == 0 ) EvalLpp();
    }
    if ( wpp ) {
      wppDeltaZ = region.z / sizeHistWpp; // computing the binSlab size (slab thickness) - used in wpp.c
      if ( ( stepCount % stepWpp ) == 0 ) EvalWpp( sizeHistWpp, wppDeltaZ );
    }
    if ( epp ) {
      eppDeltaZ = region.z / sizeHistEpp; // computing the binSlab size (slab thickness) - used in epp.c
      if ( ( stepCount % stepEpp ) == 0 ) EvalEpp();
    }
    if ( zConstraint ) {
      ComputeAndUnwrapSoluteLatCom();
      if ( stepCount == 0 || stepCount % 50 == 0 ) {
	PrintSoluteLatMotion(); 
      }
    }
    if ( stepCount >= stepLimit ) moreCycles = 0;
  } /* end  while ( moreCycles ) */  
  if ( doCheckpoint ) PutCheckpoint(); /* put final checkpoint */

  fPtr = fopen( "data.final", "w" ); 
  WriteLammpsData( fPtr );
  fclose( fPtr );

  return 0; /* make the compiler happy */
}

